Importing necessary libraries

library(tidyverse)
library(zoo)
library(xts)
library(quantmod)
library(imputeTS)
library(tseries)
library(rugarch)
library(FinTS)
library(car)
library(fBasics)

Lets load additional function defined which will be used to to easily compare information criteria for ARCH and GARCH models.

source("/Users/elgun/Desktop/ESTIMATING-VALUE-AT-RISK-OF-A-PORTFOLIO-WITH-GARCH-FAMILY-MODELS-Time-Series-project/functions/compare_ICs_ugarchfit.R")

Import the data

# Fetch data for Dow Jones
Dow_Jones <- getSymbols("^DJI", src = "yahoo", from = "2019-01-01", to = "2024-06-01", auto.assign = FALSE)
# Rename the column
colnames(Dow_Jones) <- c("Open", "High", "Low", "Close", "Volume", "Adjusted")

# Fetch data for TSMC
Tmsc <- getSymbols("TSM", src = "yahoo", from = "2019-01-01", to = "2024-06-01", auto.assign = FALSE)
# Rename the column
colnames(Tmsc) <- c("Open", "High", "Low", "Close", "Volume", "Adjusted")

# Fetch data for Pound
Pound <- getSymbols("GBPUSD=X", src = "yahoo", from = "2019-01-01", to = "2024-06-01", auto.assign = FALSE)
# Rename the column
colnames(Pound) <- c("Open", "High", "Low", "Close", "Volume","Adjusted")

# Fetch data for Silver
Silver <- getSymbols("SI=F", src = "yahoo", from = "2019-01-01", to = "2024-06-01", auto.assign = FALSE)
# Rename the column
colnames(Silver) <- c("Open", "High", "Low", "Close",  "Volume","Adjusted")

# Fetch data for BNB
BNB <- getSymbols("BNB-USD", src = "yahoo", from = "2019-01-01", to = "2024-06-01", auto.assign = FALSE)
# Rename the column
colnames(BNB) <- c("Open", "High", "Low", "Close", "Volume", "Adjusted")

print(Pound)

Data analysis and stylized facts

Handling with NA values use na_locf

Dow_Jones$Adjusted <- na_locf(Dow_Jones$Adjusted)
Tmsc$Adjusted <- na_locf(Tmsc$Adjusted)
Pound$Adjusted <- na_locf(Pound$Adjusted)
Silver$Adjusted <- na_locf(Silver$Adjusted)
BNB$Adjusted <- na_locf(BNB$Adjusted)

Log return of ecah asset

Dow_Jones$r <- diff.xts(log(Dow_Jones$Adjusted))
Tmsc$r <- diff.xts(log(Tmsc$Adjusted))
Pound$r <- diff.xts(log(Pound$Adjusted))
Silver$r <- diff.xts(log(Silver$Adjusted))
BNB$r <- diff.xts(log(BNB$Adjusted))
print(Dow_Jones$r)

Portfolio construction

# Calculate the combined returns of all assets
combined_returns <- Dow_Jones$r + Tmsc$r + Pound$r + Silver$r + BNB$r
# Calculate the portfolio returns with equal weights
weight <- 0.2
portfolio_returns <- combined_returns * weight
# Convert the xts object to a data frame
portfolio_df1 <- data.frame(date = index(combined_returns), portfolio_returns = coredata(portfolio_returns))
portfolio_df1 <- na.omit(portfolio_df1)
# Ensure the date column is in Date format
portfolio_df1$date <- as.Date(portfolio_df1$date)
# Removing observations after `2023-01-01`
portfolio_df <- portfolio_df1 %>%
  filter(date <= as.Date("2023-01-01"))

# Print the portfolio data frame
print(na.omit(portfolio_df))

Log-returns plot

plot(portfolio_df$r, 
     col = "red",
     major.ticks = "years", 
     grid.ticks.on = "years",
     main = "Log-returns of Portfolio")

Plot of ACF for returns

acf(portfolio_df$r,lag.max = 36, na.action = na.pass,col="darkblue",lwd=7,main="ACF of Portfolio returns")

acf(portfolio_df$r, lag.max = 36, na.action = na.pass, ylim = c(-0.2, 0.2),
    col = "darkblue", lwd = 7, main = "ACF of Portfolio returns")

Plot of ACF for squared returns.

acf(portfolio_df$r^2,lag.max = 100,na.action = na.pass,col="darkblue",lwd=7,main="ACF of Portfolio squared returns")

acf(portfolio_df$r^2,lag.max = 100,na.action = na.pass,
    ylim= c(0, 0.5),col="darkblue",lwd=7,main="CF of Portfolio squared returns")

Lag 1,2,3,4,6,8 seems significant.

Normality test

basicStats((portfolio_df$r))

Skewness is negative and we also observe strong excess kurtosis.

tibble(r = as.numeric(portfolio_df$r)) %>%
  ggplot(aes(r)) +
  geom_histogram(aes(y =..density..),
                 colour = "black", 
                 fill = "pink") +
  stat_function(fun = dnorm, 
                args = list(mean = mean(portfolio_df$r), 
                            sd = sd(portfolio_df$r))) +
  theme_bw() + 
  labs(
    title = "Density of the portfolio log-returns", 
    y = "", x = "",
    caption = "source: own calculations"
  )

As we can see, the distribution of the returns is highly leptokurtic.

jarque.bera.test(portfolio_df$r)

The null hypothesis about normality strongly rejected!

Testing for ARCH effects.

ArchTest(portfolio_df$r,lags = 5)

There is very strong evidence to reject the null-hypothesis that there are no ARCH effects in the data.

Build GARCH models

Estimating the GARCH(1,1)

Here, first we have to define a model specification:

spec <- ugarchspec(variance.model = list(model="sGARCH",garchOrder= c(1,1)),mean.model = list(armaOrder = c(0,0), include.mean= T),distribution.model = "norm")

Estimate the model:

portfolio_df.garch11 <- ugarchfit(spec = spec, 
                           data = portfolio_df$r)

Now, we can see the results:

portfolio_df.garch11

Interpration of results

1.Ljung-Box Tests for Residuals and Squared Residuals: The Ljung-Box tests for both standardized residuals and squared standardized residuals indicate that there is no significant autocorrelation remaining in the residuals at the tested lags.

2.Arch Test for Residuals: No remaining ARCH effects, suggesting the model has adequately captured conditional heteroscedasticity.

3.Nyblom Stability Test: Indicates potential instability in the model parameters over time, particularly the constant term (ω).

4.Sign Bias Test: Shows no evidence of asymmetry in the impact of positive and negative shocks on volatility, suggesting the model handles these shocks symmetrically.

5.Adjusted Pearson Goodness-of-Fit Test: Reveals significant deviations from normality in the model’s residuals, indicating the model may not adequately capture the true distribution of the data.

Conditional SD (vs |returns|)

plot(portfolio_df.garch11, which=3)

There are periods where the volatility is relatively low and stable, suggesting less market turbulence.

ACF of Standardized Residuals

plot(portfolio_df.garch11, which=10)

ACF of Squared Standardized Residuals

plot(portfolio_df.garch11, which=11)

Plot of ACFs support our test results.

plot(portfolio_df.garch11, which=12)

Positive and negative shocks of equal magnitude have the same impact on volatility. The U-shaped curve indicates that both positive and negative shocks increase volatility, but the impact diminishes as the magnitude of the shock decreases.

The EGARCH model

Here, first we have to define a model specification:

spec <- ugarchspec(# variance equation
                   variance.model = list(model = "eGARCH", 
                                         garchOrder = c(1, 1)),
                   # sGARCH would stand for standard GARCH model
                   # mean equation
                   mean.model = list(armaOrder = c(0, 0), 
                                     include.mean = TRUE), 
                   # assumed distribution of errors
                   distribution.model = "norm")

Estimate the model:

portfolio_df.egarch11 <- ugarchfit(spec = spec, 
                           data = portfolio_df$r)

Now, we can see the results:

portfolio_df.egarch11

Interpration

1.Ljung-Box Tests for Residuals and Squared Residuals: The Ljung-Box test results indicate no significant serial correlation, suggesting that the model is appropriate for capturing the autocorrelation in the time series data.Similarly, the test results for the squared residuals indicate no significant serial correlation, suggesting that the model adequately captures the conditional heteroscedasticity (volatility clustering).

2.Arch Test for Residuals: The fitted model adequately captures both autocorrelation and conditional heteroscedasticity in the time series data, indicating no significant serial correlation or ARCH effects in the residuals.

3.Nyblom Stability Test: The Nyblom stability test indicates that all parameters in the model are stable, with no significant evidence of instability in either individual or joint tests.

4.Sign Bias Test: Based on the Sign Bias Test results, there is no significant indication of a systematic bias in the signs of the residuals, suggesting that the model’s residuals exhibit no discernible pattern in their signs.

5.Adjusted Pearson Goodness-of-Fit Test: The Adjusted Pearson Goodness-of-Fit Test results suggest that the model does not adequately fit the observed data across all groups. The low p-values indicate significant evidence against the model’s fit, implying that the model may need adjustments or improvements to better capture the characteristics of the data.

Conditional SD (vs |returns|)

plot(portfolio_df.egarch11, which = 3)

There are periods where the volatility is relatively low and stable, suggesting less market turbulence.

ACF of Standardized Residuals

plot(portfolio_df.egarch11, which = 10)

ACF of Squared Standardized Residuals

plot(portfolio_df.egarch11, which = 11)

Plot of ACFs support our test results.

News-impact curve

plot(portfolio_df.egarch11, which = 12)

Positive and negative shocks of equal magnitude have the same impact on volatility. The U-shaped curve indicates that both positive and negative shocks increase volatility, but the impact diminishes as the magnitude of the shock decreases.

The GARCH-in-Mean model

Let’s first define a model specification:

spec <- ugarchspec(# variance equation
                   variance.model = list(model = "sGARCH", 
                                         # sGARCH = standard GARCH
                                         garchOrder = c(1, 1)),
                   # mean equation - lets turn on the intercept term
                   mean.model = list(armaOrder = c(0, 0), 
                                     include.mean = TRUE,
                       # we add an element to the mean equation,
                       # which can be either stdev (archpow 1)
                       # or var (archpow=2)
                                     archm = TRUE, archpow = 1), 
                   # assumed distribution of errors
                   distribution.model = "norm")

Then, we can estimate the model:

portfolio_df.garchm11 <- ugarchfit(spec = spec, 
                           data = portfolio_df$r)

Let’s examine the results:

portfolio_df.garchm11

Interpration

1.Ljung-Box Tests for Residuals and Squared Residuals: Both the Weighted Ljung-Box tests on standardized residuals and standardized squared residuals suggest that there is no evidence of serial correlation in the residuals at the 5% significance level.

2.Arch Test for Residuals: These results suggest that there is no significant evidence of autoregressive conditional heteroskedasticity (ARCH effects) in the residuals of the model at the 5% significance level.

3.Nyblom Stability Test: The Nyblom stability test suggests that the model’s parameters are stable.

4.Sign Bias Test: The Sign Bias Test suggests that there’s no significant systematic bias in the signs of the residuals.

5.Adjusted Pearson Goodness-of-Fit Test: Reject the null hypothesis and conclude that there is a significant difference between the observed and expected frequencies in each group.

Conditional SD (vs |returns|)

plot(portfolio_df.garchm11, which = 3)

There are periods where the volatility is relatively low and stable, suggesting less market turbulence.

ACF of Standardized Residuals

plot(portfolio_df.garchm11, which = 10)

ACF of Squared Standardized Residuals

plot(portfolio_df.garchm11, which = 11)

Plot of ACFs support our test results.

News-impact curve

plot(portfolio_df.garchm11, which = 12)

Positive and negative shocks of equal magnitude have the same impact on volatility. The U-shaped curve indicates that both positive and negative shocks increase volatility, but the impact diminishes as the magnitude of the shock decreases.

The GARCH-t model

Let’s see whether conditional distribution of the error term can be better described by the t-Student distribution.

Let’s first define a model specification:

spec <- ugarchspec(# variance equation
                   variance.model = list(model = "sGARCH", 
                                         garchOrder = c(1, 1)),
                   # mean equation
                   mean.model = list(armaOrder = c(0, 0), 
                                     include.mean = TRUE), 
                   # assumed distribution of errors
                   distribution.model = "std") # std = t-Student

Then, we estimate the model:

portfolio_df.garcht11 <- ugarchfit(spec = spec, 
                           data = portfolio_df$r)

Model summary:

portfolio_df.garcht11

Interpration

1.Ljung-Box Tests for Residuals and Squared Residuals: There is no evidence of significant serial correlation in the residuals of your time series model.

2.Arch Test for Residuals: There is no significant ARCH effect present in the residuals of your time series model at the conventional significance level.

3.Nyblom Stability Test: There is evidence of instability in the parameters of your time series model, particularly for “omega” and “beta1”.

4.Sign Bias Test: There is no evidence of systematic sign bias in the residuals of your time series model.

5.Adjusted Pearson Goodness-of-Fit Test: There isn’t significant evidence to suggest a lack of fit between the observed and expected frequencies in each group, based on the adjusted Pearson goodness-of-fit test.

Conditional SD (vs |returns|)

plot(portfolio_df.garcht11, which = 3)

There are periods where the volatility is relatively low and stable, suggesting less market turbulence.

ACF of Standardized Residuals

plot(portfolio_df.garcht11, which = 10)

ACF of Squared Standardized Residuals

plot(portfolio_df.garcht11, which = 11)

Plot of ACFs support our test results.

News-impact curve

plot(portfolio_df.garcht11, which = 12)

Positive and negative shocks of equal magnitude have the same impact on volatility. The U-shaped curve indicates that both positive and negative shocks increase volatility, but the impact diminishes as the magnitude of the shock decreases.

The EGARCH in mean-t model

Let’s first define the model specification:

spec <- ugarchspec(# variance equation
                   variance.model = list(model = "eGARCH", 
                                         garchOrder = c(1, 1)),
                   # mean equation
                   mean.model = list(armaOrder = c(0, 0), 
                                     include.mean = TRUE,
                                     archm = TRUE, archpow = 1), 
                   # assumed distribution of errors
                   distribution.model = "std") # std = t-Student

Now, we estimate the model:

portfolio_df.egarchmt11 <- ugarchfit(spec = spec, 
                             data = portfolio_df$r)

The model summary:

portfolio_df.egarchmt11

Interpration

1.Ljung-Box Tests for Residuals and Squared Residuals: There is no significant serial correlation present in the residuals or their squares based on these tests.

2.Arch Test for Residuals: There is not enough evidence to reject the null hypothesis of no conditional heteroskedasticity at any of the tested lag lengths. Therefore, based on these results, the data do not exhibit significant conditional heteroskedasticity.

3.Nyblom Stability Test: There is no evidence to suggest parameter instability in model.

4.Sign Bias Test: There’s no significant evidence to suggest a systematic bias in the predictions of model.

5.Adjusted Pearson Goodness-of-Fit Test: There’s no strong evidence to suggest significant deviations between the observed data and the specified theoretical distribution for most groups.

Conditional SD (vs |returns|)

plot(portfolio_df.garcht11, which = 3)

There are periods where the volatility is relatively low and stable, suggesting less market turbulence.

ACF of Standardized Residuals

plot(portfolio_df.garcht11, which = 10)

ACF of Squared Standardized Residuals

plot(portfolio_df.garcht11, which = 11)

Plot of ACFs support our test results.

News-impact curve

plot(portfolio_df.garcht11, which = 12)

Positive and negative shocks of equal magnitude have the same impact on volatility. The U-shaped curve indicates that both positive and negative shocks increase volatility, but the impact diminishes as the magnitude of the shock decreases.

Now, let’s compare information criteria for all models:

compare_ICs_ugarchfit(c("portfolio_df.garch11",
                        "portfolio_df.egarch11",
                        "portfolio_df.garchm11", 
                        "portfolio_df.garcht11", 
                        "portfolio_df.egarchmt11"))

Calculation of VaR

VaR in the IN-SAMPLE period for the GARCH-t model

Standardization of returns.

portfolio_df$rstd <- (portfolio_df$r - mean(portfolio_df$r, na.rm = T)) /
    sd(portfolio_df$r ,na.rm = T)
tail(portfolio_df)

1% empirical quantile

q01 <- quantile(portfolio_df$rstd, 0.01, na.rm = T)
q01

For comparison: 1% quantile of standard normal distribution

qnorm(0.01, 0, 1)
str(portfolio_df.garcht11)
head(portfolio_df.garcht11@fit$sigma)

Calculating value-at-risk (VaR) the GARCH-t model.

portfolio_df$var <- q01*portfolio_df.garcht11@fit$sigma
tail(portfolio_df)

Plot of returns vs value-at-risk.

plot(portfolio_df$date,portfolio_df$r,  col="red", 
     lwd=1,type="l", ylim= c(-0.1, 0.1))
abline(h = 0, lty = 2)
lines(portfolio_df$date,portfolio_df$var,col="green", type="l")

Post-2020, the portfolio returns (red line) show continued volatility but at a relatively stable level compared to the spike in early 2020. The VaR line (green line) gradually recovers from its early 2020 dip but remains at a lower level compared to pre-2020, indicating a sustained higher risk. The plot provides a comprehensive view of how the portfolio returns and the associated risks have evolved over time, particularly highlighting the impact of the COVID-19 pandemic on market volatility and risk levels.

In how many days losses were higher than the assumed value-at-risk?

sum(portfolio_df$r < portfolio_df$var)/ length(portfolio_df$var)

VaR in the OUT-OF-SAMPLE period fot the GARCH-t model

Plot of conditional standard deviation and its on-day ahead prediction.

plot(ugarchforecast(portfolio_df.garcht11, na.ahead=1), which=3)

The forecast indicates that volatility is expected to rise slightly in the near future.

Plot of conditional standard deviation forecasts in the long run.

plot(ugarchforecast(portfolio_df.garcht11, n.ahead = 200), which=3)

The unconditional sigma will decrease and stabilize over time.

We can combine them with the in-sample estimation of conditional standard deviation

sigma.forecast.longrun <- ugarchforecast(portfolio_df.garcht11, n.ahead = 500)
unconditional_sigma <- 
  sqrt(
    portfolio_df.garch11@model$pars["omega", 1] / 
      (1 - 
         portfolio_df.garcht11@model$pars["alpha1", 1] -
         portfolio_df.garcht11@model$pars["beta1", 1]))
plot(
  c(as.numeric(portfolio_df.garcht11@fit$sigma),
    as.numeric(sigma.forecast.longrun@forecast$sigmaFor)),
  type = "l",
  ylab = "sigma")
abline(h = unconditional_sigma, col = "red")

Yet a better idea is to annualize all of these values.

plot(
  c(as.numeric(portfolio_df.garcht11@fit$sigma * sqrt(252)),
    as.numeric(sigma.forecast.longrun@forecast$sigmaFor * sqrt(252))),
  type = "l",
  ylab = "sigma annualized")
abline(h = unconditional_sigma * sqrt(252), col = "red")

The large spike in the early part of the plot indicates a period of extreme market stress, which is common during significant financial events such as the onset of the COVID-19 pandemic.

Now, the loop below calculates predictions of VaR for the whole OUT-OF-SAMPLE period

portfolio_df1 <-
  portfolio_df1 %>%
  mutate(obs = row_number())
start   <- portfolio_df1$obs[portfolio_df1$date == as.Date("2023-01-03")]
finish  <- portfolio_df1$obs[portfolio_df1$date == as.Date("2024-05-30")]
portfolio_df2 <- portfolio_df1[start:finish, ]
VaR <- rep(NA, times = finish - start + 1)
mu     <- rep(NA, times = finish - start + 1)
omega  <- rep(NA, times = finish - start + 1)
alpha1 <- rep(NA, times = finish - start + 1)
beta1  <- rep(NA, times = finish - start + 1)

Calculation lasts for ~90-120 seconds:

time1 <- Sys.time()
for (k in start:finish) {
    tmp.data <- portfolio_df1[portfolio_df1$obs <= (k - 1), ]
    tmp.data <- tmp.data[as.Date("2019-01-01") <= tmp.data$date, ]
    tmp.data$rstd <- 
      (tmp.data$r - mean(tmp.data$r, na.rm = T)) / sd(tmp.data$r, na.rm = T)
    q01  <- quantile(tmp.data$rstd, 0.01, na.rm = T)
    spec <- 
      ugarchspec(variance.model = list(model = "sGARCH", garchOrder = c(1, 1)),
                 mean.model     = list(armaOrder = c(0, 0), include.mean = T),
                 distribution.model = "std")
    tmp.garch11           <- ugarchfit(spec = spec, data = na.omit(tmp.data$r))
    sigma.forecast        <- ugarchforecast(tmp.garch11, n.ahead = 1)
    sigma.forecast2       <- sigma.forecast@forecast$sigmaFor[1, 1]
    VaR[k - start + 1]    <- q01 * sigma.forecast2
    mu[k - start + 1]     <- tmp.garch11@fit$coef[1]
    omega[k - start + 1]  <- tmp.garch11@fit$coef[2]
    alpha1[k - start + 1] <- tmp.garch11@fit$coef[3]
    beta1[k - start + 1]  <- tmp.garch11@fit$coef[4]
  }
time2 <- Sys.time()
time2 - time1

Adding VaR values to the data.frame objects.

portfolio_df2$Var <- VaR 
head(portfolio_df2$Var)

Plot of returns vs. VaR in the OUT-OF-SAMPLE period

plot(portfolio_df2$date, portfolio_df2$r, col = "red", lwd = 1, type = 'l',
     ylim = c(-0.20, 0.20))
abline(h = 0, lty = 2)
lines(portfolio_df2$date, portfolio_df2$Var, type = 'l', col = "green")

The actual returns are more volatile than the predicted VaR. The red line frequently stays above the green line, it suggests the portfolio is performing better than the risk estimate.

In how many days losses were higher than the assumed VaR?

sum(portfolio_df2$r < portfolio_df2$Var) / length(portfolio_df2$Var)

VaR in the IN-SAMPLE period for the GARCH-mean-t model

str(portfolio_df.egarchmt11)
head(portfolio_df.egarchmt11@fit$sigma)

Calculating value-at-risk (VaR) the GARCH-mean-t model.

portfolio_df$varmt <- q01*portfolio_df.egarchmt11@fit$sigma
tail(portfolio_df)

Plot of returns vs value-at-risk.

plot(portfolio_df$date,portfolio_df$r,  col="red", 
     lwd=1,type="l", ylim= c(-0.1, 0.1))
abline(h = 0, lty = 2)
lines(portfolio_df$date,portfolio_df$varmt,col="green", type="l")

Post-2020, the portfolio returns (red line) show continued volatility but at a relatively stable level compared to the spike in early 2020. The VaR line (green line) gradually recovers from its early 2020 dip but remains at a lower level compared to pre-2020, indicating a sustained higher risk. The plot provides a comprehensive view of how the portfolio returns and the associated risks have evolved over time, particularly highlighting the impact of the COVID-19 pandemic on market volatility and risk levels.

In how many days losses were higher than the assumed value-at-risk?

sum(portfolio_df$r < portfolio_df$varmt)/ length(portfolio_df$varmt)

VaR in the OUT-OF-SAMPLE period for the GARCH-mean-t model

Plot of conditional standard deviation and its on-day ahead prediction.

plot(ugarchforecast(portfolio_df.egarchmt11, na.ahead=1), which=3)

The forecast indicates that volatility is expected to rise slightly in the near future.

Plot of conditional standard deviation forecasts in the long run.

plot(ugarchforecast(portfolio_df.egarchmt11, n.ahead = 200), which=3)

the unconditional sigma will decrease and stabilize over time.

We can combine them with the in-sample estimation of conditional standard deviation

sigma.forecast.longrun <- ugarchforecast(portfolio_df.egarchmt11, n.ahead = 500)
unconditional_sigma <- 
  sqrt(
    portfolio_df.egarchmt11@model$pars["omega", 1] / 
      (1 - 
         portfolio_df.egarchmt11@model$pars["alpha1", 1] -
         portfolio_df.egarchmt11@model$pars["beta1", 1]))
plot(
  c(as.numeric(portfolio_df.egarchmt11@fit$sigma),
    as.numeric(sigma.forecast.longrun@forecast$sigmaFor)),
  type = "l",
  ylab = "sigma")
abline(h = unconditional_sigma, col = "red")

Yet a better idea is to annualize all of these values.

plot(
  c(as.numeric(portfolio_df.egarchmt11@fit$sigma * sqrt(252)),
    as.numeric(sigma.forecast.longrun@forecast$sigmaFor * sqrt(252))),
  type = "l",
  ylab = "sigma annualized")
abline(h = unconditional_sigma * sqrt(252), col = "red")

The large spike in the early part of the plot indicates a period of extreme market stress, which is common during significant financial events such as the onset of the COVID-19 pandemic.

Now, the loop below calculates predictions of VaR for the whole OUT-OF-SAMPLE period

portfolio_df1 <-
  portfolio_df1 %>%
  mutate(obs = row_number())
start   <- portfolio_df1$obs[portfolio_df1$date == as.Date("2023-01-03")]
finish  <- portfolio_df1$obs[portfolio_df1$date == as.Date("2024-05-30")]
portfolio_df3 <- portfolio_df1[start:finish, ]
VaR <- rep(NA, times = finish - start + 1)
mu     <- rep(NA, times = finish - start + 1)
omega  <- rep(NA, times = finish - start + 1)
alpha1 <- rep(NA, times = finish - start + 1)
beta1  <- rep(NA, times = finish - start + 1)

Calculation lasts for ~90-120 seconds:

time1 <- Sys.time()
for (k in start:finish) {
    tmp.data <- portfolio_df1[portfolio_df1$obs <= (k - 1), ]
    tmp.data <- tmp.data[as.Date("2019-01-01") <= tmp.data$date, ]
    tmp.data$rstd <- 
      (tmp.data$r - mean(tmp.data$r, na.rm = T)) / sd(tmp.data$r, na.rm = T)
    q01  <- quantile(tmp.data$rstd, 0.01, na.rm = T)
    spec <- 
      ugarchspec(variance.model = list(model = "eGARCH", garchOrder = c(1, 1)),
                 mean.model     = list(armaOrder = c(0, 0), include.mean = T),
                 distribution.model = "std")
    tmp.garch11           <- ugarchfit(spec = spec, data = na.omit(tmp.data$r))
    sigma.forecast        <- ugarchforecast(tmp.garch11, n.ahead = 1)
    sigma.forecast2       <- sigma.forecast@forecast$sigmaFor[1, 1]
    VaR[k - start + 1]    <- q01 * sigma.forecast2
    mu[k - start + 1]     <- tmp.garch11@fit$coef[1]
    omega[k - start + 1]  <- tmp.garch11@fit$coef[2]
    alpha1[k - start + 1] <- tmp.garch11@fit$coef[3]
    beta1[k - start + 1]  <- tmp.garch11@fit$coef[4]
  }
time2 <- Sys.time()
time2 - time1

Adding VaR values to the data.frame objects.

portfolio_df3$Var <- VaR
head(portfolio_df3$Var)

Plot of returns vs. VaR in the OUT-OF-SAMPLE period

plot(portfolio_df3$date, portfolio_df2$r, col = "red", lwd = 1, type = 'l',
     ylim = c(-0.20, 0.20))
abline(h = 0, lty = 2)
lines(portfolio_df2$date, portfolio_df2$Var, type = 'l', col = "green")

The actual returns are more volatile than the predicted VaR. The red line frequently stays above the green line, it suggests the portfolio is performing better than the risk estimate.

In how many days losses were higher than the assumed VaR?

sum(portfolio_df3$r < portfolio_df3$Var) / length(portfolio_df3$Var)

Conclusion

The performance of both models are similar.